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We present a quantum algorithm which simulates the quantum kicked rotator model exponentially 
faster than classical algorithms. This shows that important physical problems of quantum chaos, 
localization and Anderson transition can be modelled efficiently on a quantum computer. We also 
show that a similar algorithm simulates efficiently classical chaos in certain area-preserving maps. 
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The massive parallelism of quantum evolution allows 
to manipulate simultaneously exponentially many states 
through entanglement (see reviews [§],0]). That opens 
new horizons for computations based on quantum me- 
chanics, as was stressed by Feynman [||. Nevertheless 
it is not so obvious if this parallelism can be used to 
speed up exponentially any given computational algo- 
rithm. There is certainly no systematic way to do this, 
and that is why so much interest has been generated by 
the Shor quantum algorithm Q which factorizes large 
numbers exponentially faster than any known classical 
algorithm. At present very few other quantum algo- 
rithms have been found in which substantial computa- 
tional gain is achieved compared to classical computing. 
Among them is the Grover algorithm || which signifi- 
cantly accelerates the problem of searching an unsorted 
database, although the gain is not exponential. In addi- 
tion quantum computers can be used as analog machines 
to simulate some many-body quantum systems which are 
hardly accessible in usual computer simulations S. In 
this way some systems such as spin lattices can be mod- 
elled very naturally for example by cold atoms in optical 
lattices H . However this kind of simulation is restricted 
to systems whose physical elements are related or simi- 
lar to the qubits (spin halves) of a quantum computer. 
At present apart from these natural examples there are 
no developed quantum algorithms which would allow to 
reach exponential gain in the computation of the quan- 
tum dynamics of physical systems. 

It is therefore desirable to find a quantum algorithm 
corresponding to a physical model with rich and com- 
plex quantum dynamics. During the last decades, it has 
been understood that generally the dynamics of classi- 
cal nonlinear systems is chaotic Bj. The corresponding 
quantum dynamics, called quantum chaos, demonstrates 
a rich and complex behaviour even for systems with only 
few degrees of freedom and rather simple Hamiltonians 
Pj . One of the cornerstone models in the study of quan- 
tum chaos is the kicked rotator. In the classical limit, 
this model reduces to an area-preserving map called the 
Chirikov standard map JTcfl which has applications in dif- 
ferent fields of physics, such as particle confinement in 



magnetic traps, beam dynamics in accelerators, comet 
trajectories and many others ||. The map depends on 
only one parameter, and depending on its value the sys- 
tem can be in the near-integrable bounded regime, with 
Kolmogorov-Arnold-Moser (KAM) curves, or in the fully 
chaotic regime with diffusive growth of momentum which 
statistically can be described by the Fokker-Planck equa- 
tion. In between these two regimes the phase space of 
the system has a complex hierarchical structure with in- 
tegrable islands surrounded by a chaotic sea at smaller 
and smaller scales. The quantum dynamics correspond- 
ing to these different regimes have been intensively stud- 
ied by different groups in the field of quantum chaos . 
Many phenomena of general importance are present in 
this model, including quantum ergodicity, spectral statis- 
tics as in random matrix theory, quantum KAM regime, 
and many others. However the most unexpected quan- 
tum effect is the phenomenon of dynamical localization, 
in which quantum interference suppresses chaotic dif- 
fusion in momentum, leading to exponentially localized 
eigenstates. This effect has close analogy with Anderson 
localization of electrons in disordered materials jl^] , and 
therefore this model enables to study also the properties 
of Anderson localization, a solid-state problem still un- 
der intense investigation nowadays. The quantum kicked 
rotator describes also the properties of microwave ioniza- 
tion of Rydberg atoms jljj . It has been realized exper- 
imentally with cold atoms, and the effects of dynamical 
localization, external noise and decoherence have been 
studied experimentally H]. 

In this paper we present a quantum algorithm which 
computes the evolution of the quantum kicked rotator ex- 
ponentially faster than any classical computation. It sim- 
ulates the kicked rotator with N levels in 0((log 2 N) 3 ) 
operations instead of 0(A^log 2 N) for the classical algo- 
rithm. 

The classical dynamics of the system is given by the 
Chirikov standard map 
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where (n, 6) is the pair of conjugated momentum (ac- 
tion) and angle variables, and the bars denote the re- 
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suiting variables after one iteration of the map. In this 
way the dynamics develops on a cylinder (periodicity in 
9) which can be also closed to form a torus of length 
TV = 2ttL/T where L is an integer. The classical dynam- 
ics depends only on one single chaos parameter K = kT, 
so that the motion is globally chaotic for K > 0.9716... 
For if > 1 the orbits spread diffusively in n with dif- 
fusion rate D = n 2 /t as k 2 /2 where t is measured in 
number of iterations (kicks) |lC|||- In the chaos regime, 
the dynamics is characterized by positive Kolmogorov- 
Sinai entropy h w \n(K/2) > 0, due to which trajectories 
diverge exponentially and roundoff errors grow exponen- 
tially with t p|. 

The quantum evolution during one period is described 
by an unitary operator acting on the wave function ip 



individual qubits and two-qubit gates from the ground 
state 1 0. . .0 >. For example, a typical initial state used 
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where h = —id/d9, K = 
cal limit corresponds to k 

K = kT = const Jilt] • T ne quantum interference leads 
to exponential localization of the eigenstates Xm{n) of 
the operator U in the momentum space n with envelope 
X m (rt) ~ exp(— | n — m\/l)/y/l where m marks also the 
center of the eigenstate and I is the localization length. In 
the regime of quantum chaos (k ^> K ^> 1) this length is 
determined by the classical diffusion rate I = D/2 sa fc 2 /4 
PI . The evolution takes place on N levels with periodic 
boundary conditions. For I 3> N the eigenstates become 
ergodic and the level spacing statistics is described by 
random matrix theory 11 1. Therefore depending on the 
parameters various regimes of quantum chaos can be in- 
vestigated in this single model. 

The evolution operator U is the product of two unitary 
operators Uk — exp(— ikcosQ) and Ut = exp(— iTh 2 /2) 
which represent respectively the effects of a kick and free 
rotation. These operators are diagonal in the 9 and n rep- 
resentations respectively. Due to that the most efficient 
way to simulate the quantum dynamics of this system 
on a classical computer is to perform forward/backward 
fast Fourier transforms (FFT) to go from one represen- 
tation to the other 16 1, doing diagonal multiplications 
by Uk and Ut between each FFT. In this way for a sys- 
tem with iV levels the FFT requires 0(N\og 2 (N)) oper- 
ations and the diagonal multiplications take O(N) oper- 
ations, so that evolution on one period is performed in 
0(Nlog 2 (N)) operations p7] . Our construction of the 
quantum algorithm keeps the global structure of this clas- 
sical algorithm, and uses quantum parallelism to speed 
up exponentially each algorithmic step. 

Step I: Preparation of the input state. We consider a 
system of n q qubits; the Hilbert space of dimension N = 
2 rLq is used to describe N momentum states (eigenstates 
of the operator Ut for < n < N — lin binary code) on 
which evolves the kicked rotator. An initial state i/j(0) = 
12n=o a n\ n > a * time t — is prepared by rotations of 



in the studies of the kicked rotator dynamics 1 1 1 , such 
as ip(0) = I-/V/2 >, requires only one individual rotation. 
We need also auxiliary registers which will be used later; 
at the moment they are all in the ground state. 

Step II: Action of free propagation operator Ut- In 
the n representation Ut is diagonal so that Ut\ti >= 
exp(— iTn 2 /2)\n >. The simultaneous multiplica- 
tion of the N coefficients can be done in n 2 gate 



operations. Indeed, if n 



YTj=Q l a j 2i , then 
n 2 = J2 3l J2 ®h a h 2h+h - Therefore exp(-iTn 2 /2) = 
ITj 1 .j 2 exp(— iToj 1 a! i j 2 2 3 ' 1+ - ? ' 2-1 ) with cx.j 1 2 = or 1. As a 
result, this step can be realized with n 2 operations of the 
two-qubit gate applied to each qubit pair (j'1,.72) which 
keeps the states 1 00 >, 1 01 >,|10 > unchanged while 
1 11 > is transformed to exp(-iT2-' 1 +-' 2 - :L )|ll > 

Step III: Change from n to 9 representation. In anal- 
ogy with the classical algorithm, we can use the quantum 
Fourier transform (QFT) (described in detail, for exam- 
ple in [Q). The QFT requires 0(n 2 ) operations with 
one-qubit rotations and two-qubits gates similar to the 
ones described above. After the QFT, we obtain the wave 
function in the 9 representation, X^o 1 ^tl^i >j where the 
9i are the binary codes of N — 2"« discretized angles. We 
note that the QFT was discussed in |l^| for simulating a 
rather specific model, the quantum baker map. 

Step IV: Construction of a supplementary register 
holding the cosines of angles. This step transforms 
Eto 1 b i\0i > |0 > into S^" 1 bi\9i > | cos0< >. The sec- 
ond register is of course present since step I, but is used 
only at that step. After these operations, it contains the 
binary codes of the N values of cos#i, correlated with 
9i in the first register. The number of qubits p in this 
second register sets the precision of the cosines at 2~ p , 
and should be equal or greater than n q . This register 
will be used in the next step to perform the kick oper- 
ator Uk- To realize this transformation, we need a few 
auxiliary registers which can be erased at the end. First 
we precompute the 2n q values cos(27r/2 : '), sin(27r/2 : '), for 
j = 1, ..,n q with precision 2~ p . This can be done quan- 
tum mechanically or classically in polynomial time by 
first computing the case of smallest angle and then us- 
ing recursive relations, doubling the angle each time. 
Also, other classical methods converging superexponen- 
tially (e. g. Newton's) can be used. We decompose 0, 
in binary code 6i = Ei=i Aj27r/2 J and use the formula 
exp^) = n^ 1 exp(z/3y27r/2^) = n^ 1 (cos(A,2^/2 J ) + 
i sin(/3y-27r/2- 7 )), with /3y = or 1, to compute cos#i and 
sin 9i in An q multiplications. This can be done in parallel 
for all N values of 9i in 0(n q p 2 ) gate operations. We need 
for that an auxiliary register with p qubits on which sin 9i 
is built. We start therefore with Eilo 1 > |1 > |0 >, 
then we perform the transformation EiLo 1 > l c > 
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\d >^ Eto 1 b i\ i > I cos(/3 tf 27r/2J)c-sin(Aj27r/2J')d > 
| sin(/3y27r/2^)c + cos(/3. iJ 27r/2^)d > for j = 1, ..,n q (ini- 
tially |c >= |1 >, |tf >= |0 >). This transformation 
needs a controlled multiplier with the qubit j as con- 
trol qubit. The quantum circuits realizing a controlled 
multiplier are described in [^0|, and require 0(p 2 ) gate 
operations for each multiplier. After n q transformations 
we obtain Y^i=o bi\Qi > \ cos Oi > \ sin 9i >. The to- 
tal number of gate operations for step IV is therefore 
0(n q p 2 ). We note again that p determines the preci- 
sion with which the cos Oi are computed. A reasonable 
regime is n q < p < 2n q which gives a total number of 
gate operations for this step of 0(n q ). Besides the three 
registers already described, some auxiliary registers are 
necessary to perform these operations. It can be done 
with 5 additional registers of size p, but it is probably 
possible to decrease this number. After that only the 
registers holding Oi and cosOi will be used, and all others 
can be erased. The auxiliary registers holding interme- 
diate values of cosines and sines can be reused after each 
multiplication by exp(i/9ij-27r/2 3 ) is performed, since they 
can be erased by using multiplication by exp(— ij3ij2Tr/2 J ) 
and subtraction. 

Step V: Action of kick operator Uk- After the previ- 
ous steps, the state of the system is ip = J2iLo > 
|cos0j >. In the angle representation, the action of Uk 
is diagonal so that Uk\0i >— exp(— ik cos Qi)\0i >. Each 
state \0i > is entangled with | cos^ > holding the binary 
code of cos6>, = X^=i7y2 _ "', with 7^ = or 1. Since 
exp(— ikcosOi) — n^ =1 exp(— ik"fij2~ J ), to perform the 
multiplication, it is therefore enough to apply to each 
qubit of the second register the one-qubit gate which 
takes |0 > to |0 > and |1 > to exp(— iA;2~ J ')|l >. Only 
p gate operations (with n q < p < 2n q ) are used for this 
transformation. After this | cos Oi > is reversibly erased. 
As a result of this step, the state of the system is now 
YsZTq 1 h exp(-ik cos Oi)\0i > |0 >. 

Step VI: Change from to n representation. This step 
is similar to step III, it performs backwards the QFT on 
the first register \0i > and returns the wave function to 
the momentum basis in 0(n 2 ) operations. This gives the 
wave function of the kicked rotator after one iteration of 
map (||) (one kick step). 

In this way one kick iteration requires 0(n q ) gate op- 
erations. Subsequent kicks are realized using steps II-III- 
IV-V-VI, since step I is done only once. As a result, a 
quantum computer can perform the kicked rotator evolu- 
tion exponentially faster than a classical computer (which 
needs 0(n q 2 n ") operations per kick). Several successive 
measurements after t iterations can give the largest prob- 
abilities \ai\ 2 in the momentum basis from which for ex- 
ample the localization length I can be extracted with only 
few measurements. Other average characteristics can be 
obtained efficiently by performing the QFT followed only 
by few measurements, giving for example the largest har- 



monics of the probability distribution. 

Some modifications in this algorithm are possible. In- 
stead of the additional register | cos 0% > one can use a 
register holding an uniform polynomial approximation 
P(0) w cos 6*. This is actually the technique used in 
classical computation, where Chebychev polynomials are 
used. The construction of the polynomial of degree p is 
done iteratively, starting from the lowest-degree monom 
and multiplying it by Oi to obtain the next one. An addi- 
tional register holds the current power of Oi between each 
step. The coefficients of the polynomial should be pre- 
computed in advance. The transition from one degree 
to another requires 0(p 2 ) gate operations for p qubits, 
so that the whole process needs 0(p 3 ) operations. For 
p ~ n q , this is similar to step IV. After computing the 
register \P(0i) > the next step is unchanged and performs 
multiplication by exp(— ikP(0)) 

The generalized models of kicked rotator, where k cos 
is replaced by another function V(0), are also of interest. 
In general it is not obvious if the register \V(0i) > in 
step IV can be computed in polynomial time, so uniform 
polynomial approximations can be the only way. A case 
of particular interest is V(0) = 2arctan(£' — 2kcos0), 
where E and k are parameters. In this case, the kicked 
rotator can be exactly mapped on a solid-state system on 
a chain with only nearest-neighbour hopping, as it had 
been shown in |l2] |. This computation can be done via 
uniform polynomial approximation. It may also be pos- 
sible that this V(0) can be computed directly in parallel 
from registers | cos > | sin > . 

The generalized kicked rotator can also model the lo- 
calization properties in higher dimensions. For exam- 
ple, if the parameter E varies from kick to kick, for 
example E = — 2fccos(wii) — 2/ccos(w2t) where oji/{2tt) 
and 102/ (2tt) are two incommensurate frequencies, then 
the model can be mapped on the three-dimensional An- 
derson model which displays a metal-insulator transi- 
tion pj|. For k < k c ss 0.46 the dynamics is localized 
in momentum space, while for k > k c it becomes dif- 
fusive and ergodic over all the available space of size 
N. Similar effects were observed recently in quasiperi- 
odically driven cold atoms |p3| . However, this latter 
model is still operating on a kicked rotator with one de- 
gree of freedom. One can consider a d-dimensional ver- 
sion of (||) with Uk — exp(— ik Ylt=i cos @r) an d Ut = 

exp(~iTJ2t=i n r( n r + J2t<r' n r')/2), which is directly 
related to the d-dimensional Anderson localization prob- 
lem. The simulation of this model on a system of size 
N = 2 nq for each momentum n r requires 0(N d log 2 V) 
operations on a classical computer since the total ba- 
sis contains N d levels. Even for d = 2 this model is 
hardly accessible for nowadays computers (24|. On the 
contrary, our quantum algorithm can be directly ex- 
tended to higher dimensions by increasing the number 
of registers by a factor d, and the number of operations 
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For d — 2 by changing Ut to 
igS ni ,n 2 ) it is possible to 



becomes only O(dn^). 

U T = exp(-iT{n 2 + n|)/2 
simulate the problem of two interacting particles in a lo- 
calized phase ]2j|. We note that usually one needs an 
exponential number of gates to simulate a quantum map 
operator, and a simulation in polynomial time is not ob- 
vious. For example for the extensively studied kicked top 
p6f the QFT cannot be used to change representations. 

It is interesting to note that on the basis of the algo- 
rithm constructed above, one can simulate also the clas- 
sical map ([!]). Indeed, the discretization of the map can 
be done in a symplectic way |27j on a phase space lattice 
of size N x N: 



Y = Y + S N (X); X = X + Y(modN) 



(3) 



with S N {X) = [NKsm(2%X/N)/(2w)], [...} being the 
integer part and X, Y integers. This map is area- 
preserving and invertiblc. An initial classical phase space 
density can therefore be modelled by a quantum state 

Sii,i 2 a iui2\ X h > \ Y h > 1° >• Tnen we use ste P IV 
to transform it in £) iliia a^^X^ > \Y i2 > {Sn^X^) > 

in 0(log 2 (iV) 3 ) operations, performed only once at the 
beginning. After that the map is reduced to simple 
additions in the first two registers and require only 
0(log 2 (A)) operations per iteration, simulating exponen- 
tially many trajectories in polynomial time. A simulation 
on a classical computer requires 0(N 2 ) operations per it- 
eration if one simulates a density distributed over 0(N 2 ) 
lattice cells. After t iterations aj give the density prob- 
ability distribution in phase space. QFT can be used to 
determine the main harmonics of the density, which can 
be measured by running the algorithm several times. 

In conclusion, the algorithm presented here shows that 
important systems displaying rich and complex proper- 
ties like classical and quantum chaos, quantum localiza- 
tion and ergodicity can be simulated exponentially faster 
on a quantum computer. This allows to study the time- 
evolution of these systems and the transition between 
classical and quantum mechanics in the limit of large 
quantum numbers (semiclassical limit). Using a few tens 
of qubits can probe this limit far beyond what is possible 
on nowadays computers. The problem of the precision 
of such quantum computation is crucial, and the effects 
of imperfections should be studied in detail for differ- 
ent physical properties of the simulated system (some of 
them can be quite sensitive |2q] ) . 

We thank G.Benenti and K.Frahm for stimulating dis- 
cussions. 
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